!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!
!
#if (defined (SYS_AIX) && !defined (VAR_INT64))
#define BIGENDIAN 1
#else
#define LITTLEENDIAN 1
#endif
C
C  /* Deck pckr8 */
      SUBROUTINE PCKR8(BUFFER,NDATA,BUFPCK,NBYTE,IPCKTABLE,LPACK)
*=====================================================================*
C
C     Purpose: compress BUFFER with NDATA real*8 numbers
C
C        BUFFER -- array with unpacked real*8 numbers
C        BUFPCK -- buffer with the compressed numbers
C        NDATA  -- length of unpacked array BUFFER in R*8 words
C        NBYTE  -- length of packed array BUFPCK in bytes 
C        IPCKTABLE -- packing table generated by INITPCKR8
C        LPACK  -- flag, if set to .false. PCKR8 does only a DCOPY
C        
C
C     Christof Haettig, August 1998
C
*---------------------------------------------------------------------*
      IMPLICIT NONE
C
#include "priunit.h"
C
      LOGICAL LOCDBG
      PARAMETER (LOCDBG = .FALSE.)
C
      LOGICAL LPACK
      INTEGER NBYTE, NDATA
      INTEGER IPCKTABLE(0:255)
C
      REAL*8 BUFFER(NDATA)
      REAL*8 R8WORD, ZERO
C
      PARAMETER (ZERO = 0.0D0)
C
      LOGICAL*1 BUFPCK(*)
C
      LOGICAL*1 L1BYTE(8), LTAB(4)
      INTEGER*4 I4WORD(2), ITAB 
C
      EQUIVALENCE(R8WORD,I4WORD(1))
      EQUIVALENCE(R8WORD,L1BYTE(1))
      EQUIVALENCE(ITAB,LTAB(1))
C
      INTEGER IB1, IW1, IW2, IT1, L, I, OFF, MBYTE
      INTEGER*4 M1, MASK1, MASK2, HBITS, ISIGN, ITMP
C
#ifdef LITTLEENDIAN
      PARAMETER (IT1 = 1) 
      PARAMETER (IB1 = 8) 
      PARAMETER (IW1 = 2, IW2 = 1) 
#elif BIGENDIAN
      PARAMETER (IT1 = 4) 
      PARAMETER (IB1 = 1) 
      PARAMETER (IW1 = 1, IW2 = 2) 
#endif
C
#if defined (VAR_GFORTRAN)
      ! special code to avoid gfortran out of range error
      PARAMETER (M1 = z'40000000')
      PARAMETER (MASK1 = -M1)
#else
      PARAMETER (MASK1 = z'C0000000')
#endif
      PARAMETER (MASK2 = z'3FFFFFFF')

C
C---------------------------------------------------------------
C     if packing disabled, copy input data to output array
C----------------------------------------------------------------
C
      IF (.NOT.LPACK) THEN
        CALL DCOPY(NDATA,BUFFER,1,BUFPCK,1)
        NBYTE = 8 * NDATA
        RETURN
      END IF
C
C----------------------------------------------------------------
C     Pack the input data and copy to output array
C----------------------------------------------------------------
C
#if !(defined (BIGENDIAN) || defined(LITTLEENDIAN))
      CALL QUIT('PCKR8 not implemented for this operating system.')
#endif
      NBYTE = 0
      ITAB  = 0
C
      DO I = 1, NDATA
C
         IF (LOCDBG) THEN
           WRITE (LUPRI,'(//A,I5,A,E30.20)') 'Element #',I,' = ',
     &           BUFFER(I)
           CALL R8BITREP(BUFFER(I))
         END IF
C
         R8WORD      = BUFFER(I)
C
         ISIGN       = IAND( I4WORD(IW1),MASK1)
         I4WORD(IW1) = ISHFT(I4WORD(IW1),4)
         I4WORD(IW1) = IAND( I4WORD(IW1),MASK2)
         I4WORD(IW1) = IOR(  I4WORD(IW1),ISIGN)
C
         HBITS       = ISHFT(I4WORD(IW2),-28)
         I4WORD(IW2) = ISHFT(I4WORD(IW2),4)
         I4WORD(IW1) = IOR(  I4WORD(IW1),HBITS)
C
         LTAB(IT1)   = L1BYTE(IB1)
C
         IF (LOCDBG) CALL I4BITREP(ITAB)
C
         IF (LOCDBG) CALL R8BITREP(R8WORD)
C
         MBYTE = IPCKTABLE(ITAB)
C
#ifdef LITTLEENDIAN
         OFF = 8 - MBYTE
         DO L = 1, MBYTE
           BUFPCK(NBYTE+L) = L1BYTE(OFF+L)
         END DO
#elif defined(BIGENDIAN)
         OFF = MBYTE + 1
         DO L = 1, MBYTE
           BUFPCK(NBYTE+L) = L1BYTE(OFF-L)
         END DO
#endif
C
         NBYTE = NBYTE + MBYTE
C
         IF (LOCDBG) THEN
           WRITE (LUPRI,'(3I5)') ITAB, NBYTE, MBYTE
         END IF
C
      END DO
C
C     RETURN
      END
*=====================================================================*
C  /* Deck unpckr8 */
      SUBROUTINE UNPCKR8(BUFFER,NDATA,BUFPCK,NBYTE,IPCKTABLE,LPACK)
*=====================================================================*
C
C     Purpose: unpack BUFPCK with NBYTE compressed data
C              containing NDATA real*8 numbers
C
C        BUFFER -- array with unpacked real*8 numbers
C        BUFPCK -- buffer with the compressed numbers
C        NDATA  -- length of unpacked array BUFUPK in R*8 words
C        NBYTE  -- length of packed array BUFPCK in bytes 
C        IPCKTABLE -- packing table generated by INITPCKR8
C        LPACK  -- flag, if set to .false. PCKR8 does only a DCOPY
C
C     note that both dimensions NBYTE and NDATA are input!
C
C     Christof Haettig, August 1998
C
*---------------------------------------------------------------------*
      IMPLICIT NONE
C
#include "priunit.h"      
C
      LOGICAL LOCDBG
      PARAMETER (LOCDBG = .FALSE.)
C
      LOGICAL LPACK
      INTEGER NBYTE, NDATA
      INTEGER IPCKTABLE(0:255)
C
      DOUBLE PRECISION BUFFER(NDATA)
      DOUBLE PRECISION R8WORD, ZERO
C
      PARAMETER (ZERO = 0.0D0)
C
      LOGICAL*1 BUFPCK(NBYTE)
C
      LOGICAL*1 L1BYTE(8), LTAB(4)
      INTEGER*4 I4WORD(2), ITAB
C
      EQUIVALENCE(R8WORD,I4WORD(1))
      EQUIVALENCE(R8WORD,L1BYTE(1))
      EQUIVALENCE(ITAB,LTAB(1))
C
      INTEGER IB1, IB4, IW1, IW2, IT1, L, OFF, KBYTE, MBYTE, MDATA
      INTEGER*4 M1, M4, MASK1, MASK3, MASK4, MASK5, HBITS, ISIGN, ITMP
C
#ifdef LITTLEENDIAN
      PARAMETER (IT1 = 1) 
      PARAMETER (IB1 = 8) 
      PARAMETER (IW1 = 2, IW2 = 1) 
#elif defined(BIGENDIAN)
      PARAMETER (IT1 = 4) 
      PARAMETER (IB1 = 1) 
      PARAMETER (IW1 = 1, IW2 = 2) 
#endif
C
#if defined (VAR_GFORTRAN)
      ! special code to avoid gfortran out of range error
      PARAMETER (M1    = z'40000000')
      PARAMETER (MASK1 = -M1)
      PARAMETER (MASK3 = z'3C000000')
      PARAMETER (M4    = z'3C000001')
      PARAMETER (MASK4 = -M4)
      PARAMETER (MASK5 = z'0000000F')
#else
      PARAMETER (MASK1 = z'C0000000')
      PARAMETER (MASK3 = z'3C000000')
      PARAMETER (MASK4 = z'C3FFFFFF')
      PARAMETER (MASK5 = z'0000000F')
#endif
C
C---------------------------------------------------------------
C     if packing disabled, copy input data to output array
C----------------------------------------------------------------
C
      IF (.NOT.LPACK) THEN
        CALL DCOPY(NDATA,BUFPCK,1,BUFFER,1)
        RETURN
      END IF
C
C---------------------------------------------------------------
C     Unpack the input data and copy to output array
C----------------------------------------------------------------
C
#if !(defined (BIGENDIAN) || defined(LITTLEENDIAN))
      CALL QUIT('PCKR8 not implemented for this operating system.')
#endif
      MBYTE  = NBYTE
      MDATA  = NDATA
      ITAB   = 0
      R8WORD = ZERO
C
      DO WHILE ( MDATA .GT. 0 )
C
         LTAB(IT1)   = BUFPCK(MBYTE)
C
         IF (LOCDBG) CALL I4BITREP(ITAB)
C
         KBYTE = IPCKTABLE(ITAB)
C
#ifdef LITTLEENDIAN
         MBYTE = MBYTE - KBYTE
         OFF = 8 - KBYTE
         DO L = 1, KBYTE
           L1BYTE(OFF+L) = BUFPCK(MBYTE+L) 
         END DO
#elif defined(BIGENDIAN)
         DO L = 1, KBYTE
           L1BYTE(L) = BUFPCK(MBYTE+1-L)
         END DO
         MBYTE = MBYTE - KBYTE
#endif
C
         IF (LOCDBG) CALL R8BITREP(R8WORD)
C
         HBITS = IAND(I4WORD(IW1),MASK5)
         HBITS = ISHFT(HBITS,28)
         I4WORD(IW2) = ISHFT(I4WORD(IW2),-4)
         I4WORD(IW2) = IOR(I4WORD(IW2),HBITS)
C
         ISIGN       = IAND(I4WORD(IW1),MASK1)
         I4WORD(IW1) = ISHFT(I4WORD(IW1),-4)
         I4WORD(IW1) = IOR( I4WORD(IW1),ISIGN)
C
         IF ( BTEST(ITAB,6) ) THEN
           I4WORD(IW1) = IAND(I4WORD(IW1),MASK4)
         ELSE 
           I4WORD(IW1) = IOR( I4WORD(IW1),MASK3)
         END IF
C
         BUFFER(MDATA)   = R8WORD
C
         IF (LOCDBG) THEN
           WRITE (LUPRI,'(//A,I5,A,E30.20)') 'Element #',MDATA,' = ',
     &           R8WORD
           WRITE (LUPRI,'(3I5)') ITAB, MBYTE+KBYTE, KBYTE
         END IF
C
         MDATA = MDATA - 1
C
      END DO
C
C     RETURN
      END
*=====================================================================*
C  /* Deck initpckr8 */
      SUBROUTINE INITPCKR8(THRESHOLD,IPCKTABLE,ENABLED)
*=====================================================================*
C
C     Purpose: initialize exponent table for packing/unpacking
C              of arrays with real*8 (double precision) numbers
C
C     Christof Haettig, August 1998
C
*---------------------------------------------------------------------*
      IMPLICIT NONE
C
#include "priunit.h"
C
      LOGICAL LOCDBG
      PARAMETER (LOCDBG = .FALSE.)
C
      INTEGER IPCKTABLE(0:255)
      LOGICAL ENABLED
C
      DOUBLE PRECISION THRESHOLD, ZERO, ERRMAX(8)
      DOUBLE PRECISION R8WORD0, R8WORD1, PCKTHRS
C
      PARAMETER (ZERO = 0.0D0)
C
      INTEGER IEXP, IWORD
C
      LOGICAL*1 L1BYTE0(8), L1BYTE1(8)
      INTEGER*4 I4WORD0(2), I4WORD1(2)
      INTEGER*4 IMASK(2,7), JMASK
C
      EQUIVALENCE(R8WORD0,I4WORD0(1))
      EQUIVALENCE(R8WORD0,L1BYTE0(1))
      EQUIVALENCE(R8WORD1,I4WORD1(1))
      EQUIVALENCE(R8WORD1,L1BYTE1(1))
C
      INTEGER IW1, IW2, L
C
#undef NOPACKING
#if LITTLEENDIAN
      PARAMETER (IW1 = 2, IW2 = 1) 
#elif defined(BIGENDIAN)
      PARAMETER (IW1 = 1, IW2 = 2) 
#else
#define NOPACKING
#endif
C
#if defined (LITTLEENDIAN) || defined(BIGENDIAN)
      PARAMETER (JMASK = z'3C000000')
C
#if defined (VAR_GFORTRAN)
      ! special code to avoid gfortran out of range error
      INTEGER*4 Ftmp, F0000000
      PARAMETER ( Ftmp     = z'10000000' )
      PARAMETER ( F0000000 = -Ftmp)
      DATA IMASK / z'00000000', z'000FF000',
     &             z'00000000', z'00000FF0',
     &               F0000000 , z'0000000F',
     &             z'0FF00000', z'00000000', 
     &             z'000FF000', z'00000000',
     &             z'00000FF0', z'00000000',
     &             z'0000000F', z'00000000'/
#else
      DATA IMASK / z'00000000', z'000FF000',
     &             z'00000000', z'00000FF0',
     &             z'F0000000', z'0000000F',
     &             z'0FF00000', z'00000000', 
     &             z'000FF000', z'00000000',
     &             z'00000FF0', z'00000000',
     &             z'0000000F', z'00000000'/
#endif
#else
# define NOPACKING
#endif
C
C---------------------------------------------------------------
C     check if packing can be used on this archicture:
C---------------------------------------------------------------
C
#if defined NOPACKING
      ENABLED = .FALSE.
      RETURN
#else
      ENABLED = .TRUE.
#endif     
C
C---------------------------------------------------------------
C     save packing accuracy on common block
C---------------------------------------------------------------
C
      PCKTHRS = THRESHOLD
C
C---------------------------------------------------------------
C     loop over all exponent values and find the maximum errors
C     introduced by truncating at a certain byte length
C----------------------------------------------------------------
C
      DO IEXP = 0, 255
C
         IF (LOCDBG) WRITE (LUPRI,'(//A,I5)') 'IEXP = ', IEXP
C
         IPCKTABLE(IEXP) = 8
C
         R8WORD0 = ZERO
C
         IF (IEXP.NE.0) THEN
C
           I4WORD0(IW1) = ISHFT(IEXP,24)
           IF (LOCDBG) CALL R8BITREP(R8WORD0)
C
           I4WORD0(IW1) = ISHFTC(I4WORD0(IW1),-4,30)
           IF (LOCDBG) CALL R8BITREP(R8WORD0)
C
           IF (.NOT. BTEST(I4WORD0(IW1),30) ) THEN
             I4WORD0(IW1) = I4WORD0(IW1) + JMASK
           END IF
C
         END IF
C
         IF (LOCDBG) THEN
           CALL R8BITREP(R8WORD0)
           WRITE (LUPRI,'(A,E20.10,/)') 'R8WORD0:',R8WORD0
         END IF
C
         R8WORD1 = R8WORD0
         IF (IEXP.EQ.0) THEN
           I4WORD1(IW1) = I4WORD1(IW1) + JMASK
         END IF

         DO L = 7, 1, -1
           I4WORD1(IW1) = I4WORD1(IW1) + IMASK(2,L)
           I4WORD1(IW2) = I4WORD1(IW2) + IMASK(1,L)
           ERRMAX(L) = DABS(R8WORD0-R8WORD1) 
           IF ( ERRMAX(L) .LT. PCKTHRS ) IPCKTABLE(IEXP) = L
         END DO
         ERRMAX(8) = ERRMAX(7) / 256.0D0
C
         IF (LOCDBG) THEN
           WRITE (LUPRI,'(/A)') ' #Bytes / Truncation error:'
           WRITE (LUPRI,'(8I10)') (L,L=1,8)
           WRITE (LUPRI,'(8E10.2)') (ERRMAX(L),L=1,8)
           WRITE (LUPRI,'(A,I5,A,8E10.2)') ' selected length:',
     &          IPCKTABLE(IEXP),
     &                    ' -- accepted error:',ERRMAX(IPCKTABLE(IEXP))
         END IF
C
      END DO
C
C     RETURN
      END
*=====================================================================*
      SUBROUTINE R8BITREP(R8WORD)
*---------------------------------------------------------------------*
C
C     Purpose: print bit representation of a real*8 number
C
C     Christof Haettig, August 1998
C
*---------------------------------------------------------------------*
      IMPLICIT NONE
C
#include "priunit.h"
C
      REAL*8    R8WORD, R8COPY
      INTEGER*8 K, J, I8COPY
      CHARACTER*64 STRING

      EQUIVALENCE (R8COPY,I8COPY)

      R8COPY = R8WORD

      DO J = 1, 64
         IF ( BTEST(I8COPY,J-1) )  THEN
            STRING(J:J) = '1'
         ELSE
            STRING(J:J) = '0'
         END IF
      END DO

      WRITE (LUPRI,'(4(1X,A16))') STRING(01:16), STRING(17:32),
     &                     STRING(33:48), STRING(49:64)

      RETURN
      END 
*=====================================================================*
      SUBROUTINE I4BITREP(I4WORD)
*---------------------------------------------------------------------*
C
C     Purpose: print bit representation of a i*4 number
C
C     Christof Haettig, August 1998
C
*---------------------------------------------------------------------*
      IMPLICIT NONE
C
#include "priunit.h"      
C
      INTEGER*4 I4WORD, I4COPY
      INTEGER J
      CHARACTER*32 STRING

      I4COPY = I4WORD

      DO J = 1, 32
        IF ( BTEST(I4COPY,J-1) )  THEN
           STRING(J:J) = '1'
        ELSE
           STRING(J:J) = '0'
        END IF
      END DO

      WRITE (LUPRI,'(2(1X,A16))') STRING(1:16), STRING(17:32)
      WRITE (LUPRI,*) I4COPY

      RETURN
      END 
*=====================================================================*
